﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace new_problem
{
    internal class DataGen
    {
        public static readonly Int32 N = 50;
        public static readonly Double sigma = -0.6;
        public static readonly Double h = 15.0;
        public static readonly Double F = 6.67384; //*Math.Pow(10, -20);

        public static Double SurfacePoint(Double x, Double y)
        {
            x -= 5;
            y -= 5;

            if (x <= 170 && x >= 80 && y <= 170 && y >= 80)
            {
                return -h + 2;
            }
            else return -h;
        }

        //начальная аппроксимация
        public static Double[,] StartApproximation()
        {
            Double[,] result = new Double[N + 1, N + 1];
            for (int i = 0; i <= N; i++)
                for (int j = 0; j <= N; j++)
                    result[i, j] = h+0.000001;
            return result;
        }

        public static Double[] GetX()
        {
            Double[] result = new Double[N + 1];
            result[0] = 5.0;
            for (int i = 1; i <= N; i++)
                result[i] = result[i - 1] + 5;
            return result;
        }
        
        public static Double[] GetY()
        {
            Double[] result = new Double[N + 1];
            result[0] = 5.0;
            for (int i = 1; i <= N; i++)
                result[i] = result[i - 1] + 5;
            return result;
        }
        
        public static Double[,] GetSurfaceGrid(Double[] x, Double[] y)
        {
            Double[,] result = new Double[N + 1,N + 1];
            for (int i = 0; i <= N; i++)
                for (int j = 0; j <= N; j++)
                    result[i, j] = SurfacePoint(x[i], y[j]);
            return result;
        }

        public static Double[,] GetFieldGrid(Double[] x, Double[] y, Double[,] surf)
        {
            double hx = (x[1] - x[0]) / N;
            double hy = (y[1] - y[0]) / N;
            double multConst = -1.0 * F * sigma * hx * hy;

            double k = 0.0;
            double[,] A = new double[N + 1, N + 1];
            //Double[,] surf = GetSurfaceGrid(x, y);

            for (int i = 0; i <= N; i++)
                for (int j = 0; j <= N; j++)
                {
                    k = IntSimpson(x[i], y[j], surf, x, y);
                    A[i, j] = multConst * k;
                }
            return A;
        }

        public static Double[,] GetFieldGrid(Double[] x, Double[] y)
        {
            double hx = (x[1] - x[0]) / N;
            double hy = (y[1] - y[0]) / N;
            double multConst = -1.0 * F * sigma * hx * hy;

            double k = 0.0;
            double[,] A = new double[N + 1, N + 1];
            Double[,] surf = GetSurfaceGrid(x, y);

            for (int i = 0; i <= N; i++)
                for (int j = 0; j <= N; j++)
                {
                    k = IntSimpson(x[i], y[j], surf, x, y);
                    A[i, j] = multConst * k;
                }
            return A;
        }

        // численное интегрирование, формула Симпсона для двойного интеграла
        private static double IntSimpson(double x1, double y1, double[,] z, double[] x, double[] y)
        {
            double f = 0;

            for (int i = 1; i < N; i += 2)
                for (int j = 1; j < N; j += 2)
                    f += Func(x[i + 1], y[j + 1], x1, y1, z[i + 1, j + 1]) + Func(x[i + 1], y[j - 1], x1, y1, z[i + 1, j - 1]) +
                        Func(x[i - 1], y[j + 1], x1, y1, z[i - 1, j + 1]) + Func(x[i - 1], y[j - 1], x1, y1, z[i - 1, j - 1]) +
                        4 * (Func(x[i], y[j + 1], x1, y1, z[i, j + 1]) + Func(x[i], y[j - 1], x1, y1, z[i, j - 1]) +
                        Func(x[i + 1], y[j], x1, y1, z[i + 1, j]) + Func(x[i - 1], y[j], x1, y1, z[i - 1, j])) +
                        16 * Func(x[i], y[j], x1, y1, z[i, j]);
            
            return f / 9.0;
        }

        // подынтегральная функция
        private static Double Func(Double x, Double y, Double x1, Double y1, Double z)
        {
            return (1 / Math.Sqrt((x - x1) * (x - x1) + (y - y1) * (y - y1) + h * h)) -
                (1 / Math.Sqrt((x - x1) * (x - x1) + (y - y1) * (y - y1) + z * z));           
        }

        // численное интегрирование, прямоугольники
        private static Double IntRect(Double[] x, Double[] y, Double x1, Double y1, Double[,] z)
        {
            double k = 0;
            for (int i = 0; i <= N; i++)
                for (int j = 0; j <= N; j++)
                    k += Func(x[i], y[j], x1, y1, z[i, j]);

            return k;
        }

        public static Double[,] GetFieldGridBySurf(Double[] x, Double[] y, Double[,] z)
        {
            double hx = (x[1] - x[0]) / N;
            double hy = (y[1] - y[0]) / N;
            double multConst = -1.0 * F * sigma * hx * hy;

            double k = 0.0;
            Double[,] A = new double[N + 1, N + 1];

            for (int i = 0; i <= N; i++)
                for (int j = 0; j <= N; j++)
                {
                    k = IntRect(x, y, x[i], y[j], z);
                    A[i, j] = multConst * k;
                }
            return A;
        }

        public static Double SumDerivative(Double[] x, Double[] y, Double x1, Double y1, Double[,] z)
        {
            double k = 0;
            for (int i = 0; i <= N; i++)
                for (int j = 0; j <= N; j++)
                    k += z[i, j]/Math.Pow((x[i] - x1)*(x[i] - x1) + (y[i] - y1)*(y[i] - y1) + z[i, j]*z[i, j], 1.5);

            return k;
        }

        public static Double[,] GetFieldDerivative(Double[] x, Double[] y, Double[,] z)
        {
            double hx = (x[1] - x[0]) / N;
            double hy = (y[1] - y[0]) / N;
            double multConst = F * sigma * hx * hy;

            double k = 0.0;
            Double[,] Ad = new double[N + 1,N + 1];

            for(int i = 0; i <= N; i++)
                for(int j = 0; j <= N; j++)
                {
                    k = SumDerivative(x, y, x[i], y[j], z);
                    Ad[i, j] = multConst*k;
                }
            return Ad;
        }

        /**
         * Стабилизирующий функционал - сумма частных производных по всем точкам сетки
         */ 
        public static Double StabilisingFunctional(Double[,] z)
        {
            Double fdx = 0;
            Double fdy = 0;
            Double[,] f = new Double[N + 3, N + 3];
            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < N; j++)
                {
                    f[i, j] = z[i, j];
                }
            }
            for (int i = 0; i < N; i++)
            {
                f[i, N + 1] = f[i,N];
                f[i, N + 2] = f[i, N];
                f[N + 2, i] = f[N, i];
                f[N + 1, i] = f[N, i];
            }

            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < N; j++)
                {
                    fdx += (-3 * f[i, j] + 4 * f[i + 1, j] - f[i + 2, j]) / 4 * h;
                    fdy += (-3 * f[i, j] + 4 * f[i, j + 1] - f[i, j + 2]) / 4 * h;
                }
            }
            return fdy + fdx;
        }
    }
}
